#include <istream>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <cstdlib>

#include <pthread.h>

#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/filesystem/path.hpp>

#include <print_error.h>
#include <tree_UPGMA.h>
#include <tree.h>
#include <tree_ops.h>
#include <portable_timer.h>
#include <progress_indicator.h>
#include <matrix.h>
#include <matrix_io.h>
#define DEBUG_LEVEL 6
#include <debug_func.h>

#include "get_arg.h"

using boost::format;
using boost::lexical_cast;
using boost::filesystem::path;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::cerr;
using std::ifstream;
using std::map;
using std::pair;
using std::istream;
using std::max;
using std::min;
using namespace bcpp_matrix;

struct t_incorrect_pairs
{
	string cluster_name;
	unsigned cnt_genes;
	unsigned cnt_missing_pairs;
	unsigned cnt_expected_pairs;
	t_incorrect_pairs(const t_matrix_f& data, const string& cluster_name_):
		cluster_name(cluster_name_),
		cnt_genes(data.size()),
		cnt_missing_pairs(data.get_count_missing()),
		cnt_expected_pairs(data.get_count_pairs())
	{}
};



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

int main (int argc, char *argv[])
{
	// we use cout to print progress so flush per call
//	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;




	unsigned cnt_cluster = 0;
	unsigned cnt_subclusters = 0;
	//
	//	save clusters with missing entries
	//
	vector<t_incorrect_pairs> problems;



	//
	//	Multiple clusters can be read in with pairwise input format
	//
	while (cin.good())
	{
		++cnt_cluster;


		//
		//	read data from input file
		//
		unsigned cnt_duplicate_pairs;
		t_matrix_f initial_matrix;
		{
			args.VERBOSE() << "\tReading gene pairs and single linkage clustering...\n";
			t_progress_indicator dots(args.VERBOSE(), 20000);
			bcpp_matrix::read_matrix(cin,  initial_matrix,  args.input_format, cnt_duplicate_pairs, dots);
			if (cnt_duplicate_pairs)
				args.VERBOSE() << "\t" << cnt_duplicate_pairs << " duplicate pairs found.\n";
		}

		//
		//	apply threshold
		//
		args.VERBOSE() << "\tDiscarding values greater than " << args.threshold << " ...\n";
		for (unsigned i = 1, size = initial_matrix.size(); i < size; ++i)
			for (unsigned j = 0; j < i; ++j)
				if (initial_matrix.matrix[i][j] > args.threshold)
				{
					initial_matrix.erase_element(i, j);
					initial_matrix.erase_element(j, i);
				}


		//
		//	single linkage cluster
		//
		vector<t_matrix_f>	data;
		initial_matrix.get_single_linkage_clusters(data);

		if (!data.size())
			break;

		args.VERBOSE() << "\tWriting matrix...\n";
		for (unsigned i = 0; i < data.size(); ++i)
		{
			++cnt_subclusters;
			t_matrix_f& sub_cluster = data[i];
			sub_cluster.matrix_name = sub_cluster.matrix_name.length() ?
										sub_cluster.matrix_name :
										("Cluster# " + lexical_cast<string>(cnt_cluster)) ;

			sub_cluster.matrix_name += "." + lexical_cast<string>(i+1);

			// save problems
			if (sub_cluster.get_count_missing())
				problems.push_back(t_incorrect_pairs(sub_cluster, sub_cluster.matrix_name));

			//
			// prepare output files with specified prefices
			//
			args.open_output_files(cnt_subclusters);

			args.VERBOSE() << "\t>" << sub_cluster.matrix_name << "\n";
			args.VERBOSE() << "\t" << sub_cluster.size() << " rows/columns in the matrix.\n";

			{
				t_progress_indicator dots(args.VERBOSE(), 1);
				dots.use_timer();
				if (!args.header)
					sub_cluster.matrix_name = "";
				args.f_matrix << t_matrix_ostrm_wrapper_f(sub_cluster, args.output_format,  dots);
				//write_matrix(args.f_matrix,  sub_cluster,  args.output_format,  false, dots);
			}

		}
		if (args.input_format != bcpp_matrix::eMIO_PAIRS)
			break;
	}


	//
	//	print clusters with missing entries
	//
	if (problems.size())
	{
		std_print_warning("Some matrixes had missing pairs:");
		for (unsigned i = 0; i < problems.size(); ++i)
		{

			args.VERBOSE() << format("\t\t%|1$4|: "
									"%|2$4| genes, "
									"%|3$4| gene pairs expected, "
									"%|4$4| gene pairs missing.\n ")
					% problems[i].cluster_name
					% problems[i].cnt_genes
					% problems[i].cnt_expected_pairs
					% problems[i].cnt_missing_pairs;
		}
	}
	args.VERBOSE() << "\n\tAll finished!!\n\n";
	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(),  e.what());
		return false;
	}

	return 0;
}

